Viking Results, 2021-04

# Check requisite packages are installed.
packages <- c(
  "plotly", 
  "dplyr"
)
for (pkg in packages) {
  library(pkg, character.only = TRUE)
}

What do the results look like?

dirViking <- file.path(
  getwd(), "LCAB_LawMorton1996-NumericalPoolCommunityScaling"
)
dirVikingResults <- file.path(
  dirViking, "results-2021-04"
)
resultFormat <- paste0(
  "run-", 
  "%d", # Combination Number, or CombnNum.
  "-", 
  "%s", # Run Seed.
  ".RDS"
)

# Copied from LawMorton1996-NumericalPoolCommunityScaling-Calculation.R
# TODO: In the future, make this a separate file that everyone can call...
set.seed(38427042)

basal <- c(3, 10, 30, 100, 300, 1000)
consumer <- c(3, 10, 30, 100, 300, 1000) * 2
events <- (max(basal) + max(consumer)) * 2
runs <- 100

logBodySize <- c(-2, -1, -1, 1) # Morton and Law 1997 version.
parameters <- c(0.01, 10, 0.5, 0.2, 100, 0.1)

# Need to rerun seedsPrep to get the random number generation right for seedsRun
seedsPrep <- runif(2 * length(basal) * length(consumer)) * 1E8
seedsRun <- runif(runs * length(basal) * length(consumer)) * 1E8

paramFrame <- with(list(
  b = rep(basal, times = length(consumer)),
  c = rep(consumer, each = length(basal)),
  s1 = seedsPrep[1:(length(basal) * length(consumer))],
  s2 = seedsPrep[
    (length(basal) * length(consumer) + 1):(
      2 * length(basal) * length(consumer))
  ],
  sR = seedsRun
), {
  temp <- data.frame(
    CombnNum = 0,
    Basals = b,
    Consumers = c,
    SeedPool = s1,
    SeedMat = s2,
    SeedRuns = "",
    SeedRunsNum = 0,
    EndStates = I(rep(list(""), length(b))),
    EndStatesNum = 0,
    EndStateSizes = I(rep(list(""), length(b))),
    EndStateAssembly = I(rep(list(""), length(b)))
  )
  for (i in 1:nrow(temp)) {
    seeds <- sR[((i - 1) * runs + 1) : (i * runs)]
    temp$SeedRuns[i] <- toString(seeds) # CSV
    temp$SeedRunsNum[i] <- length(seeds)
  }
  temp$CombnNum <- 1:nrow(temp)
  temp
})

# Note: n + 2 end states. Failure to finish, failure to obtain state, and state.
for (i in 1:nrow(paramFrame)) {
  resultsList <- list(
    "No Run" = 0,
    "No State" = 0
  )
  resultsSize <- list(
    "0" = 0
  )
  resultsAssembly <- list(
    "No Run" = data.frame(),
    "No State" = data.frame()
  )
  seeds <- unlist(strsplit(paramFrame$SeedRuns[i], ', '))
  for (seed in seeds) {
    fileName <- file.path(
      dirVikingResults,
      sprintf(resultFormat, paramFrame$CombnNum[i], seed)
    )
    
    if (file.exists(fileName)) {
      temp <- load(fileName)
      temp <- eval(parse(text = temp)) # Get objects.
      
      if (is.data.frame(temp)) {
        community <- toString(
          temp[[ncol(temp)]][[nrow(temp)]]
        )
        size <- toString(length(temp[[ncol(temp)]][[nrow(temp)]]))
        
        if (community == "") {
          resultsList$`No State` <- resultsList$`No State` + 1
          resultsSize$`0` <- resultsSize$`0` + 1
          
        } else if (community %in% names(resultsList)) {
          resultsList[[community]] <- resultsList[[community]] + 1
          resultsSize[[size]] <- resultsSize[[size]] + 1
          
        } else {
          resultsList[[community]] <- 1
          resultsAssembly[[community]] <- temp
          
          if (size %in% resultsSize) {
            resultsSize[[size]] <- resultsSize[[size]] + 1
          } else {
            resultsSize[[size]] <- 1
          }
        }
      } else {
        resultsList$`No State` <- resultsList$`No State` + 1
        resultsSize$`0` <- resultsSize$`0` + 1
      }
    } else {
      resultsList$`No Run` <- resultsList$`No Run` + 1
      resultsSize$`0` <- resultsSize$`0` + 1
    }
  }
  
  paramFrame$EndStates[[i]] <- resultsList
  paramFrame$EndStatesNum[i] <- length(resultsList) - 2 # ! No State, No Run
  paramFrame$EndStateSizes[[i]] <- resultsSize
  paramFrame$EndStateSizesNum[i] <- length(resultsSize) - 1 # ! 0
  paramFrame$EndStateAssembly[[i]] <- resultsAssembly
}

# X, Y, Basal and Consumer.
# Z = Sizes of the Endstates.

plotScalingData <- data.frame(
  CombnNum = rep(paramFrame$CombnNum, paramFrame$EndStatesNum),
  Basals = rep(paramFrame$Basals, paramFrame$EndStatesNum),
  Consumers = rep(paramFrame$Consumers, paramFrame$EndStatesNum)
)

# Communities
comms <- unlist(lapply(paramFrame$EndStates, names))
freqs <- unlist(paramFrame$EndStates)
asmbl <- unlist(paramFrame$EndStateAssembly, recursive = FALSE)
asmbl <- asmbl[comms != "No Run" & comms != "No State"]
freqs <- freqs[comms != "No Run" & comms != "No State"]
comms <- comms[comms != "No Run" & comms != "No State"]

asmbl <- lapply(asmbl, function(d) {
  d %>% dplyr::filter(Result.Outcome != "Type 1 (Failure)" & 
                        Result.Outcome != "Present")
})

plotScalingData$Communities <- comms
plotScalingData$CommunityFreq <- freqs
plotScalingData$CommunitySeq <- asmbl

# Community Size
temp <- unlist(lapply(strsplit(plotScalingData$Communities, ','), length))
plotScalingData$CommunitySize <- temp

# For usage by the reader.

plotScaling <- plotly::plot_ly(
  plotScalingData,
  x = ~Basals,
  y = ~Consumers,
  z = ~CommunitySize
)

plotScaling <- plotly::add_markers(plotScaling)

plotScaling <- plotly::layout(
  plotScaling,
  scene = list(
    xaxis = list(type = "log"),
    yaxis = list(type = "log")
  )
)

plotScaling

plotScalingData %>% dplyr::select(-CommunitySeq)

How do they compare to each other?

# > runif(1) * 1E8
# [1] 82598679
set.seed(82598679)

temp <- load(file.path(
  dirViking, 
  "LawMorton1996-NumericalPoolCommunityScaling-PoolMats.RDS"
))
mats  <- eval(parse(text = temp[1]))
pools <- eval(parse(text = temp[2]))

candidateData <- plotScalingData %>% dplyr::group_by(
  CombnNum
) %>% dplyr::mutate(
  OtherSteadyStates = dplyr::n() - 1
) %>% dplyr::filter(
  OtherSteadyStates > 0
)
candidateData %>% dplyr::select(-CommunitySeq)

candidateData$CommunityAbund <- ""

Warning messages:
1: Unknown or uninitialised column: `CommunityProd`. 
2: Unknown or uninitialised column: `CommunityProd`. 
3: Unknown or uninitialised column: `CommunityProd`. 

for (r in 1:nrow(candidateData)) {
  candidateData$CommunityAbund[r] <- toString(with(
    candidateData[r, ], {
      invasions <- CommunitySeq[[1]]$Result.Addition
      abundance <- rep(0, Basals + Consumers)
      for (i in invasions) {
        abundance[i] <- abundance[i] + 1
        abundance <- RMTRCode2::quiet(
          deSolve::ode(
            #rootSolve::steady(
            y = abundance,
            times = c(0:10000),
            func = RMTRCode2::GeneralisedLotkaVolterra,
            parms = list(a = mats[[CombnNum]],
                         r = pools[[CombnNum]]$ReproductionRate),
            #method = "runsteady"
          )[10001, -1])
      }
      #print(paste("Expected", Communities))
      #print(paste("Calculated", toString(which(abundance > 0))))
      #print(paste("> Epsilon:", toString(which(abundance > 1E-6))))
      if (Communities == toString(which(abundance > 1E-6))) {
        abundance[abundance > 1E-6]
      } else {
        "Failure"
      }
    }
  ))
}

candidateData$CommunityProd <- NA
for (r in 1:nrow(candidateData)) {
  candidateData$CommunityProd[r] <- with(candidateData[r, ], 
    RMTRCode2::Productivity(
      Pool = pools[[CombnNum]], 
      InteractionMatrix = mats[[CombnNum]], 
      Community = Communities, 
      Populations = CommunityAbund
    )
  )
}

islandFUN <- function(i, dat, pool, mat, dmat) {
      temp <- dat[i, ]
      RMTRCode2::IslandDynamics(
        Pool = pool,
        InteractionMatrix = mat,
        Communities = c(
          list(temp$Communities[1]),
          rep("", nrow(dmat) - 2),
          temp$Communities[2]
        ),
        Populations = c(
          list(temp$CommunityAbund[1]),
          rep("", nrow(dmat) - 2),
          list(temp$CommunityAbund[2])
        ),
        DispersalPool = 0.0001,
        DispersalIsland = dmat,
      )
    }

# For each group,
# For each pair,
# Run Island Dynamics,
# Save the result with its pairing

for (grp in unique(candidateData$CombnNum)) {
  candidateDataSubset <- candidateData %>% dplyr::filter(CombnNum == grp)
  pairingResults<- combn(
    nrow(candidateDataSubset), 2, 
    islandFUN,
    dat = candidateDataSubset, 
    pool = pools[[grp]],
    mat = mats[[grp]],
    dmat = matrix(c(0, 1, 1, 0), nrow = 2, ncol = 2),
    simplify = FALSE
  )
}

for (grp in unique(candidateData$CombnNum)) {
  candidateDataSubset <- candidateData %>% dplyr::filter(CombnNum == grp)
  pairingResults<- combn(
    nrow(candidateDataSubset), 2, 
    islandFUN,
    dat = candidateDataSubset, 
    pool = pools[[grp]],
    mat = mats[[grp]],
    dmat = matrix(c(
      0, 1, 0, # Island 2 -> 1
      1, 0, 1, # Island 1 -> 2, Island 3 -> 2
      0, 1, 0  # Island 2 -> 3
    ), nrow = 3, ncol = 3, byrow = TRUE),
    simplify = FALSE
  )
}